c $Id: tabinit.f,v 1.14 2003/05/02 13:14:46 weber Exp $
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine loadwtab (io)
c
c     Revision : 1.0
c
coutput   : information in common-block comwid.f
c
c     load the tabulated branching ratios from disk
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      implicit none

      include 'comres.f'
      include 'comwid.f'

      integer ios, nsp, io, ver!, i
c      character*35 pwdcmd
      character*10 deftab
      character*8 defexe
      logical b

c set the defaultname of the file, containing the table
      parameter (deftab='tables.dat', defexe='uqmd.exe')

      b=io.eq.1

c get the name of the table from the environment variable
      call getenv('URQMD13_TAB',tabname)
c if it is empty, use the default name
      if (tabname(1:4).eq.'    ') then
         tabname=deftab
      endif

      if(b)write (6,*) 'Looking for the tabulated decay width...'
     &                 ,tabname
c open the table
      open (unit=75,iostat=ios,file=tabname,form='unformatted',
     .      status='old')
c if it fails ...
      if (ios.ne.0) then
c close the file handle
c         close (unit=75, status='delete')
         if(b)write (6,*) 'No file:',tabname,'in this directory'
c get the full path of the executable, ...
c         call getenv('_',pwdcmd)
c         write (6,*) 'pwd:',pwdcmd
         write (6,*) 'tabname:',tabname
c extract the path
c         i=max(index(pwdcmd,defexe),2)
c         write (tabname,*) pwdcmd (1:i-1),deftab
c         tabname=tabname(2:)
         if(b)write (6,*) 'Looking for ',tabname,'...'
c and look for a table in the directory of the executable
         open (unit=75,iostat=ios,file=tabname,
     .         form='unformatted',status='old')
      endif
c if the last 'open' command succeeds read the file
      if (ios.eq.0) then
         if(b)write (6,*) 'O.K.'
         if(b)write (6,*) 'reading...'
c read all tables
         read (unit=75, iostat=ios) ver, nsp, tabx, fbtaby, pbtaby,
     .          fmtaby, pmtaby, bwbarnorm, bwmesnorm,
     .            tabxnd, frrtaby
c caution! the file is unformatted, therefor it is system dependent!
         if(b)write (6,*) 'version=',ver
c if no errors occur ...
         if (ios.eq.0) then
            if(b)write (6,*) 'O.K.'
            wtabflg=3
c check, if the version number is correct
            if (ver.eq.tabver) then
               if(b)write (6,*) 'tabver=',ver,'  O.K.'
            else
               write (6,*) 'wrong table!'
               write (6,*) 'tabver should be',tabver,',instead of',ver
               wtabflg=0
            endif
c check, if the table has the correct 'widnsp'
            if (nsp.eq.widnsp) then
               if(b)write (6,*) 'widnsp=',nsp,'  O.K.'
            else
               write (6,*) 'wrong table!'
               write (6,*) 'widnsp should be',widnsp,', instead of',nsp
                     wtabflg=0
               endif
c if table is O.K. close file
            if (wtabflg.eq.3) then
               close (unit=75, status='keep')
c otherwise ...
            else
c delete the present table
               close (unit=75, status='delete')
               tabname=deftab
c and calculate a new one
               call mkwtab
            endif
c in case of read errors ...
         else
c delete the present table
            close (unit=75, status='delete')
            write (6,*) 'Error while reading ',tabname
            tabname=deftab
c and calculate a new one
            call mkwtab
         endif
c in any other case ...
      else
ctp         tabname=deftab
c calculate an new table
         call mkwtab
      endif

      return
      end





cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine savewtab
c
c     Revision : 1.0
c
c     save the tabulated branching ratios to disk
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      implicit none

      include 'comres.f'
      include 'comwid.f'

      integer ios

      write (6,*) 'Writing new table...'

c try to generate a new file
      open (unit=75,iostat=ios,file=tabname,form='unformatted',
     .      status='new')
c if it succedds ...
      if (ios.eq.0) then
c write the tables into the file
         write (unit=75, iostat=ios) tabver, widnsp, tabx, fbtaby,
     .        pbtaby, fmtaby, pmtaby, bwbarnorm, bwmesnorm,
     .            tabxnd, frrtaby
         if (ios.eq.0) write (6,*) 'O.K.'
c otherwise complain
      else
         write (6,*) 'Error: ',tabname,'exists!'
      endif
c close the file
      close (unit=75, status='keep')

      return
      end





cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine mkwtab
c
c     Revision : 1.0
c
coutput   : information in common-block comwid.f
c
c     tabulate the mass dependent branching ratios
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      implicit none

      include 'comres.f'
      include 'comwid.f'

      real*8 fwidth,m,first,last,delta,abl0,abln,mir,mminit,fbrancx
      real*8 massit,bran,smass,bwnorm,fppfit
      integer i,bchan,itp,isoit,cmin,cmax,i1,i2,i3,i4,ii1

      write (6,*) 'Generating table...'
c this indicates, that all tables are still empty
      wtabflg=0

c high precision splines from mintab to maxtab1
c lower precicision between maxtab1 and maxtab2

c now fill the x-values
c start with 'mintab'
      first=mintab
c 66 % of all fixpoints between mintab and maxtab1
c calculate the steps
      delta=(maxtab1-mintab)/((widnsp-1d0)*2d0/3d0)
      if (delta.le.0d0) then
         write(*,*)'(E) Please allow maxtab1>mintab in comwid'
         stop
      endif
c store the values into 'tabx'
      do 10 i=1,int(widnsp*2./3.)
         m=first+(i-1)*delta
         tabx(i)=m
 10   continue
c 33 % of all fixpoints with larger delta between maxtab1 and maxtab2
      delta=(maxtab2-maxtab1)/((widnsp-1d0)*1d0/3d0)
      if (delta.le.0d0) then
         write(*,*)'(E) Please allow maxtab2>maxtab1 in comwid'
         stop
      endif
c store the values into 'tabx'
        do 11 i=int(widnsp*2./3.)+1,widnsp
         m=maxtab1+(i-1-int(widnsp*2./3.))*delta
         tabx(i)=m
 11   continue

c now fill the y-values of the full branching ratios

c these are the first derivatives at the first an the last point
c of the interpolating function. a value greater than 1E30 signals the
c 'spline' routine to set the boundary condition for a natural spline
c with zero second derivative
      abl0=2D30
      abln=2D30

c loop over all baryons
      do 20 itp=minbar,maxbar
c loop over all x-values
         do 21 i=1,widnsp
c store the values ...
            fbtaby (i,itp,1)=fwidth(itp,isoit(itp),tabx(i))
 21      continue
c calculate the second derivate and store it in 'fbtaby(,,2)'
         call spline (tabx(1),fbtaby(1,itp,1),widnsp,abl0,abln,
     .                fbtaby(1,itp,2))
 20   continue
      write (6,*) '(1/7) ready.'

c loop over all mesons
      do 30 itp=minmes,maxmes
c loop over all x-values
         do 31 i=1,widnsp
c store the values ...
            fmtaby (i,itp,1)=fwidth(itp,isoit(itp),tabx(i))
 31      continue
c calculate the second derivate and store it in 'fmtaby(,,2)'
         call spline (tabx(1),fmtaby(1,itp,1),widnsp,abl0,abln,
     .                fmtaby(1,itp,2))
 30   continue
      write (6,*) '(2/7) ready.'

c the flag indicates, that now all full widths are tabulated
      wtabflg=1

c now fill the y-values of the partial branching ratios

c loop over all baryons
      do 40 itp=minbar,maxbar
c get the mass of this particle
         mir=massit(itp)
c get the range of possible decay channels
         call brange (itp, cmin, cmax)
c check, if there are any decay channels
         if (cmax.gt.0) then
c loop over all decay channels
            do 41 bchan=cmin,cmax
c now get the outgoing particles 'i1' and 'i2' for the channel 'j'
c 'bran' is the mass independent branching ratio (tabulated in blockres)
c 'bflag' indicates, if 'i1', 'i2' or both are broad
               call b3type (itp,bchan,bran,i1,i2,i3,i4)
c check, if decay is allowed

               smass=mminit(i2)
               if(i3.ne.0) smass=smass+mminit(i3)
               if(i4.ne.0) smass=smass+mminit(i4)

               if (bran.gt.1d-9.and.mir.gt.mminit(i1)+smass) then
c loop over all x-values
                  do 42 i=1,widnsp
c store the values

cdh
*        write(*,*)'mkwtab: i,itp,bchan=',i,itp,bchan
*       write(*,*)'mkwtab: isoit,tabx,bran,i1,i2,i3,i4=',
*    &    isoit(itp),tabx(i),bran,i1,i2,i3,i4
cdh

                     pbtaby(i,1,itp,bchan)=
     .                    fbrancx (bchan,itp,isoit(itp),tabx(i),
     .                    bran,i1,i2,i3,i4)
 42               continue
c calculate the second derivate and store it in 'pbtaby(,2,,)'
                  call spline (tabx(1),pbtaby(1,1,itp,bchan),widnsp,
     .                         abl0,abln,pbtaby(1,2,itp,bchan))
               end if
 41         continue
         end if
 40   continue
      write (6,*) '(3/7) ready.'

c loop over all mesons
      do 50 itp=minmes,maxmes
c get the mass of this particle
         mir=massit(itp)
c get the range of possible decay channels
         call brange (itp, cmin, cmax)
c check, if there are any decay channels
         if (cmax.gt.0) then
            do 51 bchan=cmin,cmax
c now get the outgoing particles 'i1' and 'i2' for the channel 'j'
c 'bran' is the mass independent branching ratio (tabulated in blockres)
c 'bflag' indicates, if 'i1', 'i2' or both are broad
               call b3type(itp,bchan,bran,i1,i2,i3,i4)
c!!!
               smass=mminit(i2)
               if(i3.ne.0) smass=smass+mminit(i3)
               if(i4.ne.0) smass=smass+mminit(i4)

               if (bran.gt.1d-9.and.mir.gt.mminit(i1)+smass) then
c loop over all x-values
                  do 52 i=1,widnsp
                     pmtaby(i,1,itp,bchan)=
     .                    fbrancx (bchan,itp,isoit(itp),tabx(i),
     .                    bran,i1,i2,i3,i4)
 52               continue
c calculate the second derivate and store it in 'pmtaby(,2,,)'
                  call spline (tabx(1),pmtaby(1,1,itp,bchan),widnsp,
     .                         abl0,abln,pmtaby(1,2,itp,bchan))
               end if
 51         continue
         end if
 50   continue

      write (6,*) '(4/7) ready.'


c calculate the norm integral of the Breit-Wigner functions
c   with mass dependent widths

c..baryons
        do 60 i=minbar,maxbar
           bwbarnorm(i)=bwnorm(i)
60        continue
      write (6,*) '(5/7) ready.'

c.. mesons
        do 61 i=minmes,maxmes
           bwmesnorm(i)=bwnorm(i)
61        continue
      write (6,*) '(6/7) ready.'

c now all branching ratios and BW-integrals are tabulated
      wtabflg=2

ce tabulate fppfit
c fill the x-values
c range of tabulated cross sections
      first=2d0*massit(nucleon)+massit(pimeson)
        last=maxtab1
c calculate the steps
c the energies are weighted quadratically
      delta=(last-first)/((widnsp-1)*2./3.)**2
c store the values into 'tabx'
c 66 % of all fixpoints between mintab and maxtab1
      do 69 i=1,int(widnsp*2./3.)
         m=first+(i-1)**2*delta
         tabxnd(i)=m
 69   continue
c 33 % of all fixpoints with larger, constant delta between maxtab1 and maxtab2
        delta=(maxtab2-last)/((widnsp-1)*1./3.)
        do 70 i=int(widnsp*2./3.)+1,widnsp
         m=maxtab1+(i-1-int(widnsp*2./3.))*delta
         tabxnd(i)=m
 70   continue


c.. all pp-exit channels
c loop over first out-particle N & D
        do 81 ii1=1,2
          if(ii1.eq.1)i1=minnuc
          if(ii1.eq.2)i1=mindel
c loop over second out-particle N(1440)..maxdel
          do 82 i2=minnuc+1,maxdel
c loop over all x-values
          do 83 i=1,widnsp
c store the values ...
              frrtaby(i,1,ii1,i2)=fppfit(99,tabxnd(i),i1,i2)
83            continue
c calculate the second derivate and store it in 'frrtaby(,,2)'
          call spline (tabxnd(1),frrtaby(1,1,ii1,i2),widnsp,abl0,abln,
     .                frrtaby(1,2,ii1,i2))
82          continue
81        continue


      write (6,*) '(7/7) ready.'

c pp cross sections are now tabulated
        wtabflg=3

c save the table on disk
      call savewtab

      return
      end



cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real*8 function splint (xa,ya,y2a,n,x)
c
c     Unit     : general infrastructure
c     Author   : (C) Copr. 1986-92 Numerical Recipes Software
c     Date     : 03/07/96
c     Revision : 1.1
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      implicit none

      include 'comres.f'
      include 'comwid.f'

      integer n
      integer k,khi,klo
      real*8 x,y,xa(n),y2a(n),ya(n)
      real*8 a,b,h

      save khi,klo
      data klo/1/
      data khi/2/

      if(khi.le.n.and.x.ge.xa(klo).and.x.lt.xa(khi))then
      elseif(khi+1.le.n.and.x.ge.xa(klo+1).and.x.lt.xa(khi+1))then
        klo=klo+1
        khi=khi+1
      else

      klo=1
      khi=n
1     if (khi-klo.gt.1) then
         k=(khi+klo)/2d0
         if(xa(k).gt.x)then
            khi=k
         else
            klo=k
         endif
      goto 1
      endif
      endif

      h=xa(khi)-xa(klo)
      if (h.eq.0.) pause 'bad xa input in splint'
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+((a*a*a-a)*y2a(klo)+
     .            (b*b*b-b)*y2a(khi))*(h*h)/6d0
      splint=y

      return
      end



cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real*8 function splintth (xa,ya,y2a,n,x,th)
c
c     Unit     : general infrastructure
c     Author   : (C) Copr. 1986-92 Numerical Recipes Software
c                modified my H. Weber
c     Date     : 03/07/96
c     Revision : 1.1
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     split routine with nice threshold behaviour for cross sections
c

      implicit none

      include 'comres.f'
      include 'comwid.f'

      integer n
      integer k,khi,klo
      real*8 x,y,xa(n),y2a(n),ya(n)
      real*8 a,b,h,th

      save khi,klo
      data klo/1/
      data khi/2/

      if(khi.le.n.and.x.ge.xa(klo).and.x.lt.xa(khi))then
      elseif(khi+1.le.n.and.x.ge.xa(klo+1).and.x.lt.xa(khi+1))then
        klo=klo+1
        khi=khi+1
      else

      klo=1
      khi=n
1     if (khi-klo.gt.1) then
         k=(khi+klo)/2d0
         if(xa(k).gt.x)then
            khi=k
         else
            klo=k
         endif
         goto 1
      endif
      endif
      h=xa(khi)-xa(klo)
      if (h.eq.0.) pause 'bad xa input in splint'
      if (xa(khi).lt.(th+2*h)) then
c linear approximation close to threshold (within 2h)
         splintth=ya(khi)*(x-th)/(xa(khi)-th)
      else
         a=(xa(khi)-x)/h
         b=(x-xa(klo))/h
         y=a*ya(klo)+b*ya(khi)+((a*a*a-a)*y2a(klo)+
     .        (b*b*b-b)*y2a(khi))*(h*h)/6d0
         splintth=y
      endif

      return
      end
